#ifndef TURBO_FUNCTIONS
#define TURBO_FUNCTIONS
//File that contains useful functions
#include <map>
#include <string>
#include <vector>
#include <fstream>
#include <sstream>
#include <iostream>
#include <boost/algorithm/string.hpp>
#include "grid.h"

#include <mkl_lapacke.h>

typedef std::pair<std::string, std::string> param_item;
typedef std::multimap<std::string, std::string> param_list;
class function
{
public:
	int size;
	std::vector<double> x, value;
	function(int _size)
	{
		size = _size;
		x.resize(size);
		value.resize(size);
	};
	function(){};
	void ReadData1(char* filename);	//first type of data arrangement in file
	void ReadData2(char* filename);	//second type of data arrangement in file
	void WriteData(char* filename);	//write our function in a file
};
void function::ReadData1(char* filename)
{
	std::ifstream ifs(filename);
	double read;
	while(!ifs.eof())
	{
		ifs >> read;
		x.push_back(read);		//coordinate first
		ifs >> read;
		value.push_back(read);
	};
	ifs.close();
	size = x.size();
	return;
};
void function::ReadData2(char* filename)
{
	std::ifstream ifs(filename);
	double read;
	while(!ifs.eof())
	{
		ifs >> read;
		value.push_back(read);		//value first
		ifs >> read;
		x.push_back(read);
	};
	size = x.size();
	ifs.close();
	return;
};
void function::WriteData(char* filename)
{
	std::ofstream ofs(filename);
	for(int i=0; i<x.size(); i++)
	{
		ofs << x[i] << " " << value[i] << "\n";
	};
	ofs.close();
	return;
};

//Output lapack matrix
void PrintBandedMatrix(int matrix_order, int n, int m, int kl, int ku) {
	if (matrix_order == LAPACK_COL_MAJOR) {
		for (int i = 0; i<n; i++) {
			for (int j = 0; j<m; j++) {
				 
				
			};
			printf("\n");
		};
	};
};

//Get parameter value from map or throw an error
std::vector<std::string> GetParam(param_list& list, std::string name) {
	std::vector<std::string> result(0);
	param_list::iterator it;	
	it = list.find(name);
	if (it != list.end()) {
		std::pair<param_list::iterator, param_list::iterator> prange = list.equal_range(name);
		for (it = prange.first; it!=prange.second; it++) result.push_back(it->second);
		return result;
	} else {
		std::cerr<<"Error: Parameter <<"<<name<<" wasnt defined.";
		exit(0);
		return result;
	};
};

//Function that parses configuration file and give the list of parameters presented there
param_list ParseConfig(const char *fname)
{
	param_list result;
	std::ifstream ifs(fname);
	std::string line;
	std::vector<std::string> strs;
	if (ifs.is_open())
	{
		while ( ifs.good() )
		{
			getline(ifs, line);
			//Parse line
			if (line.size() == 0) continue;	//If its empty line
			boost::split(strs, line, boost::is_any_of("#"));	//Get rid of comments
			line = strs[0];
			boost::trim(line);
			if (line.size() == 0) continue;	//If its empty line
			boost::split(strs, line, boost::is_any_of("="));
			std::string name = strs[0];
			std::string value = strs[1];
			boost::trim(name);
			boost::trim(value);
			/*
			if (value.size() == 0) {
				std::cerr << "Wrong syntaxis in config file "<<fname<<"\n";
				exit(0);
			};
			*/
			//Add pair to result
			//result[name] = value;
			result.insert(param_item(name, value));
		};
		ifs.close();
	} else {
		std::cerr << "Unable to open configuration file "<<fname<<"\n";
		exit(0);
	};
	return result;
};


Grid GenGrid(int N)
{
	Grid g;
	g.nodes.resize(N);
	for (int i = 0; i<N; i++) {
		Node new_node;
		new_node.GlobalIndex = i;
		new_node.P.x = (1.0*i) / (N-1);
		new_node.P.y = 0;
		new_node.P.z = 0;
		g.nodes[i] = new_node; 
	};
	g.cells.resize(N-1);	
	for (int i = 0; i<(N-1); i++) {
		Cell c;
		c.GlobalIndex = i;
		c.Nodes.resize(2);
		c.Nodes[0] = i;
		c.Nodes[1] = i+1;
		c.CellCenter = 0.5 * (g.nodes[i].P+g.nodes[i+1].P);
		c.CellVolume = (g.nodes[i].P-g.nodes[i+1].P).mod();
		c.Faces.resize(2);
		c.Faces[0] = i;
		c.Faces[1] = i+1;
		g.cells[i] = c;
	};
	g.faces.resize(N);
	for (int i = 0; i<N; i++) {				
		Face new_face;		
		new_face.GlobalIndex = i;
		if (i == N-1) {
			new_face.FaceCell_1 = i-1;
		} else {
			new_face.FaceCell_1 = i;
		};
		if ((i != 0) && (i!= (N-1))) {
			new_face.isExternal = false;
			new_face.FaceCell_2 = i-1;
		} else {
			new_face.isExternal = true;
			if (i == 0) new_face.BCMarker = 1;
			if (i == N-1) new_face.BCMarker = 2;
			new_face.FaceCell_2 = -1;
		};				
		new_face.FaceNodes.clear();
		new_face.FaceNodes.push_back(i);		
		//Calc surface normal and center position given all nodes
		new_face.FaceCenter = g.nodes[new_face.FaceNodes[0]].P;		
		new_face.FaceSquare = 1;
		if (i != N-1) {
			new_face.FaceNormal = Vector(-1,0,0);
		} else {
			new_face.FaceNormal = Vector(1,0,0);
		}
		g.faces[i] = new_face;
	};
	g.refresh_indexes();
	return g;
};


//Compress 
Grid GenGrid(int N, int M, double size_x, double size_y, double q_x, double q_y)
{
	Grid g;
	N++;
	M++;
	std::vector<double> x_p(N);
	std::vector<double> y_p(M);
	double h_x = 1.0 * size_x * (1.0 -q_x) / (1.0 - pow(q_x, N-1));
	double h_y = 1.0 * size_y * (1.0 -q_y) / (1.0 - pow(q_y, M-1));
	for (int i = 0; i<N; i++) {
		if (q_x == 1.0) {
			x_p[i] = (1.0*size_x*i) / (N-1);
		} else {
			x_p[i] = h_x * (1.0 - pow(q_x, i)) / (1.0-q_x);
		};
	};
	x_p[0] = 0;
	x_p[N-1] = size_x;
	for (int i = 0; i<M; i++) {
		if (q_y == 1.0) {
			y_p[i] = (1.0 *size_y*i) / (M-1);
		} else {
			y_p[i] = 1.0 * h_y * (1.0 - pow(q_y, i)) / (1.0-q_y);
		};
	};
	y_p[0] = 0;
	y_p[M-1] = size_y;
	g.nodes.resize(N*M);
	for (int i = 0; i<N; i++) {
		for (int j = 0; j<M; j++) {
			Node new_node;
			new_node.GlobalIndex = i + j * N;
			new_node.P.x = x_p[i];//(size_x*i) / (N-1);
			new_node.P.y = y_p[j];//(size_y*j) / (M-1);
			new_node.P.z = 0;
			g.nodes[new_node.GlobalIndex] = new_node; 
		};
	};
	g.cells.resize((N-1)*(M-1));	
	for (int i = 0; i<(N-1); i++) {
		for (int j = 0; j<(M-1); j++) {
			Cell c;
			c.GlobalIndex = i + j*(N-1);
			c.Nodes.resize(4);
			c.Nodes[0] = i + j*N;
			c.Nodes[1] = i + (j+1)*N;
			c.Nodes[2] = (i+1) + (j+1)*N;
			c.Nodes[3] = (i+1) + j*N;			
			c.CellCenter = Vector(0,0,0);
			for (int k = 0; k<c.Nodes.size(); k++) {
				c.CellCenter =  c.CellCenter + g.nodes[c.Nodes[k]].P;
			};
			c.CellCenter = (1.0/c.Nodes.size()) * c.CellCenter;		
			c.CellVolume = ((g.nodes[c.Nodes[1]].P-g.nodes[c.Nodes[0]].P).mod()) * ((g.nodes[c.Nodes[2]].P-g.nodes[c.Nodes[1]].P).mod());
			c.Faces.resize(4);
			c.Faces[0] = i*(M-1) + j;
			c.Faces[1] = (i+1)*(M-1) + j;
			c.Faces[2] = (i + j*(N-1)) + N*(M-1);
			c.Faces[3] = (i + (j+1)*(N-1)) + N*(M-1);			
			c.CellHSize = 0;
			g.cells[c.GlobalIndex] = c;
		};
	};
	g.faces.resize(N*(M-1) + M*(N-1));
	//Vertical first
	for (int i = 0; i<N; i++) {
		for (int j = 0; j<M-1; j++) {
			Face new_face;		
			new_face.GlobalIndex = i*(M-1) + j;
			new_face.FaceNodes.clear();
			new_face.FaceNodes.push_back(i + j * N);
			new_face.FaceNodes.push_back(i + (j+1) * N);
			new_face.FaceSquare = (g.nodes[new_face.FaceNodes[1]].P - g.nodes[new_face.FaceNodes[0]].P).mod();
			new_face.FaceCenter = Vector(0,0,0);
			for (int k = 0; k<new_face.FaceNodes.size(); k++) {
				new_face.FaceCenter = new_face.FaceCenter + g.nodes[new_face.FaceNodes[k]].P;
			};
			new_face.FaceCenter = (1.0/new_face.FaceNodes.size()) * new_face.FaceCenter;
			if (i == 0) {
				//Left border
				new_face.FaceNormal = Vector(-1, 0, 0);
				new_face.isExternal = true;
				new_face.BCMarker = 1;
				new_face.FaceCell_1 = i + j*(N-1);
				new_face.FaceCell_2 = -1;
			};
			if (i == N-1) {
				//Right border
				new_face.FaceNormal = Vector(1, 0, 0);
				new_face.isExternal = true;
				new_face.BCMarker = 2;
				new_face.FaceCell_1 = (i-1) + j*(N-1);
				new_face.FaceCell_2 = -1;
			};
			if ((i!=0) && (i!=N-1)) {
				//Inner face
				new_face.FaceNormal = Vector(1, 0, 0);
				new_face.isExternal = false;
				new_face.BCMarker = -1;
				new_face.FaceCell_1 = (i-1) + j*(N-1);
				new_face.FaceCell_2 = i + j*(N-1);
			};
			new_face.FaceNormal = new_face.FaceNormal * new_face.FaceSquare;
			g.faces[new_face.GlobalIndex] = new_face;
		};
	};
	//Horisontal now
	for (int i = 0; i<N-1; i++) {
		for (int j = 0; j<M; j++) {
			Face new_face;		
			new_face.GlobalIndex = N*(M-1) + j * (N-1) + i;
			new_face.FaceNodes.clear();
			new_face.FaceNodes.push_back(i + j * N);
			new_face.FaceNodes.push_back(i+1 + j * N);
			new_face.FaceSquare = (g.nodes[new_face.FaceNodes[1]].P - g.nodes[new_face.FaceNodes[0]].P).mod();
			new_face.FaceCenter = Vector(0,0,0);
			for (int k = 0; k<new_face.FaceNodes.size(); k++) {
				new_face.FaceCenter = new_face.FaceCenter + g.nodes[new_face.FaceNodes[k]].P;
			};
			new_face.FaceCenter = (1.0/new_face.FaceNodes.size()) * new_face.FaceCenter;
			if (j == 0) {
				//Bottom border
				new_face.FaceNormal = Vector(0, -1, 0);
				new_face.isExternal = true;
				new_face.BCMarker = 3;
				new_face.FaceCell_1 = i + j*(N-1);
				new_face.FaceCell_2 = -1;
			};
			if (j == M-1) {
				//Top border
				new_face.FaceNormal = Vector(0, 1, 0);
				new_face.isExternal = true;
				new_face.BCMarker = 4;
				new_face.FaceCell_1 = i + (j-1)*(N-1);
				new_face.FaceCell_2 = -1;
			};
			if ((j!=0) && (j!=M-1)) {
				//Inner face
				new_face.FaceNormal = Vector(0, 1, 0);
				new_face.isExternal = false;
				new_face.BCMarker = -1;
				new_face.FaceCell_1 = i + (j-1)*(N-1);
				new_face.FaceCell_2 = i + j*(N-1);
			};
			new_face.FaceNormal = new_face.FaceNormal * new_face.FaceSquare;
			g.faces[new_face.GlobalIndex] = new_face;
		};
	};	
	for (int i = 0; i<g.nodes.size(); i++) g.nodes[i].GlobalIndex++;
	for (int i = 0; i<g.cells.size(); i++) 
		for (int j = 0; j<g.cells[i].Nodes.size(); j++) g.cells[i].Nodes[j]++;
	for (int i = 0; i<g.faces.size(); i++) 
		for (int j = 0; j<g.faces[i].FaceNodes.size(); j++) g.faces[i].FaceNodes[j]++;
	g.refresh_indexes();
	return g;
};


//General exception
class Exception {
	std::string msg;
public:
	Exception(std::string m) {
		msg = m;
	};

	std::string ShowMessage() {
		std::cerr << msg.c_str() << "\n";
		return msg;
	};
};

//Calculate mach number from conservative variables values
double Mach(std::vector<double> U, double gamma) {
	if (U.size() != 5) throw new Exception("Cant calculate mach number.");
	Vector v(U[1]/U[0], U[2]/U[0], U[3]/U[0]);	
	double P = (gamma - 1.0) * (U[4] - (U[1]*U[1]+U[2]*U[2]+U[3]*U[3])/(2.0*U[0]));
	double ro = U[0];
	double c2 = (gamma * P / ro);
	//double c = sqrt(c2);
	//double c02 = cd2 + (gamma - 1) * vd.mod() * vd.mod()/2.0;
	return sqrt(v*v / c2);
};

/* Auxiliary routine: printing a matrix */
void print_matrix( char* desc, int m, int n, double* a, int lda ) {
        int i, j;
        printf( "\n %s\n", desc );
        for( i = 0; i < m; i++ ) {
                for( j = 0; j < n; j++ ) printf( " %6.6lg", a[i+j*lda] );
                printf( "\n" );
        }
};

inline double max(double a, double b) {
	return (a>b)?a:b;
};

inline double min(double a, double b) {
	return (a<b)?a:b;
};

inline int max(int a, int b) {
	return (a>b)?a:b;
};

inline int min(int a, int b) {
	return (a<b)?a:b;
};

#endif